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ABSTRACT 

Context. Solar eruptions and high flare activity often accompany the rapid rotation of sunspots. The study of sunspot rotation and the 
mechanisms driving this motion are therefore key to our understanding of how the solar atmosphere attains the conditions necessary 
for large energy release. 

Aims. We aim to demonstrate and investigate the rotation of sunspots in a 3D numerical experiment of the emergence of a magnetic 
flux tube as it rises through the solar interior and emerges into the atmosphere. Furthermore, we seek to show that the sub-photospheric 
twist stored in the interior is injected into the solar atmosphere by means of a deflnitive rotation of the sunspots. 

Methods. A numerical experiment is performed to solve the 3D resistive magnetohydrodynamic (MHD) equations using a 
Lagrangian-Remap code. We track the emergence of a toroidal flux tube as it rises through the solar interior and emerges into the 
atmosphere investigating various quantities related to both the magnetic held and plasma. 

Results. Through detailed analysis of the numerical experiment, we And clear evidence that the photospheric footprints or sunspots 
of the flux tube undergo a rotation. Signiflcant vertical vortical motions are found to develop within the two polarity sources after 
the held emerges. These rotational motions are found to leave the interior portion of the held untwisted and twist up the atmospheric 
portion of the held. This is shown by our analysis of the relative magnetic helicity as a signiflcant portion of the interior helicity is 
transported to the atmosphere. In addition, there is a substantial transport of magnetic energy to the atmosphere. Rotation angles are 
also calculated by tracing selected fleldlines; the fleldlines threading through the sunspot are found to rotate through angles of up to 
353° over the course of the experiment. We explain the rotation by an unbalanced torque produced by the magnetic tension force, 
rather than an apparent effect. 
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1. Introduction 


Sunspot rotation has attracted the interest of many researchers 
over the years from both an observational and modelling per¬ 
spective. The wide scope in observations of sunspot rotation 
merits a study of the mechanisms driving this motion. Coronal 
mass ejections (CMEs) and solar flares are often related to the 
rapid rotation of sunspots. Hence, the study of sunspot rotation 
is crucial to our understanding of such explosive events on the 
Sun as a possible mechanism for allowing the corona to achieve 
the necessary conditions for high energy release. 


Just over a century ago, whilst working at the Kodaikanal So¬ 
lar Observatory, |Evershed| ( | 1 909| ) first discovered evidence of the 
rotation of sunspots based on spectral observations. Since this 
initial evidence, sunspot rotations have been analysed in numer¬ 
ous obs ervational s t udies. In recent years, de taile d case studies 
such a s |Yan & Qu| ( |2Q07| ), |Yan et al.] (|2009|), and |Min & Chae 


P009| ) and statistical analyses from |Brown et al.|p003| ) and |Yan 
|et al. ( |2008| ) have investigated this phenomena. These studies, 
along with several others, have shown that sunspots can exhibit 
significant rotation of the order of several hundreds of degrees 
over a few days. |Yan et al.| ( [2008] ) conducted a statistical study of 
rotating sunspots using TRACE (Transition Region And Coronal 
Explorer), Hinode, and MDI (Michelson Doppler Imager) mag¬ 
netograms from 1996 to 2007, in which they individually anal¬ 
ysed 2959 active regions. They found 182 significantly rotating 
sunspots within 153 active regions. This is equivalent to approx¬ 
imately 5% of active regions harbouring rotating sunspots. On 


the other hand, [Brown et al.| ( |2003| ) conducted a more detailed 
study of seven sunspots using white light TRACE data to find 
sunspot centres and track notable features over time to calculate 
rotation rates. Rotation angles of between 40° and 200° were ob¬ 
served over periods of three to five days, resulting in an average 
rotation rate of a few degrees per hour. Six of these rotating spots 
resulted in subsequent flaring activity and the energisation of the 
corona. 


Interestingly, measurements of sunspot rotation have given 
variable results depending on the methodology employed. Eor 
example, [Min & Chaej p009| ) and [Yan et al^ p009| ) analysed 
the same a ctive region, NQAA 10930, and found notably differ¬ 
ent results. [Min & Cha^ ( |2009j ) noted a counter-clockwise rota¬ 
tion of 540 degrees over five days, whereas [Yan et alT] ( |2009[ ) 
focused on the X3.2 flare that followed the rapi d rotation of 
259 degrees over four days. Eurthermore, although [Brown et aT 


( 2003| ) and [Yan & Qu[ ( [2007) ) both conclude d that different parts 
of sunspots often rotate at differing speeds, [Brown et aL] ( [2003[ ) 
noted that the highest rotation rate was found in the penumbra 
while [Yan & Qu[p007[ ) concluded that the highest rotation rate 
was found in the umbra. This suggests that sunsp ots do not nec¬ 
essarily rotate as a rigid body. [Yan & Qu| ( [2007 ) concluded that 
twist can be created by a variation in rotation rate with distance 
from the centre of a sunspot. This twist can then be injected into 
the corona kinking the magnetic loops and driving flare activity. 

Eruptions and flares are, in fact, often correlated with ro¬ 
tational motions of sunspots. Observations have shown direct 
evidence of the energisation of the corona by these rotations. 


Article number, page 1 of 14 





























































A&A proofs: manuscript no. Sturrock2015 


Th e initiation of CME s by sunspot rotation was studied in detail 
by |Torok et al.| ( [2Q14] ) from both an observational and modelling 
viewpoint. They found that the rotation of sunspots can signifi¬ 
cantly weaken the magnetic tension of the overlying field in ac¬ 
tive regions and can then trigger an eruption in this region. This 
is an alternative explanation for eruptions caused by rotation of 
sunspots to the common theory that erupt ions are triggered by a 
direct injection of twist ( |Yan & Qu|2007| ). 

This paper aims to discuss and simulate a possible mech¬ 
anism for the rotation of sunspots. Two mechanisms were dis¬ 
cussed in [Brown et al.| ( [2003] ), namely, photospheric flows and 
magnetic flux emergence. Photospheric flows are primarily due 
to the large scale effect of differential rotation and localised mo¬ 
tions resulting from magneto-convective dynamics. The effects 
of differential rotation are kept to a minimum in [Brown et al. 


( 2003[ ) as the images were derotated prior to measuring the ve¬ 
locities. The second mechanism discussed was the emergence 
of magnetic flux. [Brown et al.[ ( [2003[ ) suggested that the pho¬ 
tospheric footpoints of a flux tube are observed to rotate as the 
tube emerges. This is the proposed mechanism for the rotation 
of sunspots which we will investigate. This mechanism seems 


the most appropriate according to th e findings of [Brown et al 
P003[ ) as well as a case study from [Min & Chae[p009 ). The 
latter study supports flux emergence as a viable mechanism as 
they discovered that the rotation speed increases in close rela¬ 
tion to the growth of the sunspot of interest, which we attribute 
to the emergence of flux. The next logical question is, how can 
flux emergence drive sunspot rotation? Two possible explana¬ 
tions were suggested by [Min & Chae[ ( |2()09[ ). They conjectured 
that the rotation may be an apparent motion caused by a twisted 
flux tube rising vertically and the fieldlines successively cross¬ 
ing the photospheric boundary at different locations as a result of 
the the twisted structure of the field. In this case there is no real 
horizontal motion of the plasma; instead this rotation is a virtual 
effect ca used by continued d isplacement of the footpoints. Alter¬ 
natively, |Min^Ch^(|20^ proposed that the observed rotation 
may represent the real horizontal motion of the plasma due to a 
net torque originating in the interior driving the plasma to rotate 
on the photospheric boundary. The torque will be examined later 
in an effort to determine its origin. 

The process of magnetic flux emergence has been consid¬ 
ered numerically in various experiments in recent years (for in¬ 
stance [FS| ( |20^ , [Archontis et ak] ( |2004| ), [HSbd et al.[ ( [20091 ) 
amongst many others). The widely accepted picture of sunspot 
formation is that an Q-shaped flux tube rises from the base of the 
convection zone until its apex intersects the photosphere to form 
a pair of sunspots. If a magnetic flux tube is in pressure balance 
and thermal equilibrium with its surroundings, the tube will be 
less dense than its surroundings, and will therefore be buoyant. 
This mechanism allows a flux tube to rise to the stably stratified 
photosphere where it remains until the tube is able to enter the 
atmosphere by initiation of a second instability, namely the mag¬ 
netic buoyancy instability. A pair of concentrations of opposite 
polarities, known as bipolar sunspots, mark the intersection of 
the field at the surface. In these experiments, the computational 
domain typically models a region from the top of the solar in¬ 
terior to the lower corona. A magnetic flux tube is then placed 
in the solar interior and is made buoyant by either introducing a 
density deficit or imposing an initial upward velocity. A recent 
review of simulations relating to this emergence process was un¬ 
dertaken by [Hood et al.[ ( [2012| ). 

Investigations of magnetic flux emergence as a possible 
mechanism for sunspot rotation have been conducted in the 
past. [Longcope & Welsc^ ( [2000| ) suggested that the rotation of 


sunspots is a consequence of the transport of helicity from the 
convection zone to the corona as a twisted flux tube emerges. 
They also demonstrated that a torsional Alfven wave will propa¬ 
gate downwards at the instance of emergence due to a mismatch 
in current between the highly twisted interior and stretched coro¬ 
nal portion of the field. In addition, [Gibson et al.[ p004[ ) ex¬ 
plained the rotation as an observational consequence of the emer¬ 
gence of a flux tube through the photosphere. An investigation 
of the transport of magnetic energy and helicity in an emerg¬ 
ing flux model was conducted by [Magara & Longcope[ ( |2003| ). 
They found that emergence generates horizontal flows as well 
as vertical flows, both of which contribute to the injection of 
energy and helicity to the atmosphere. The contributions by ver¬ 
tical flows are dominant initially but horizontal flows are the pri¬ 
mary source at a later stage. A more comprehensive study of 
the horizontal flows at the photosphere during emergence was 
later performed by [Magara[ ( [2006[ ) in which they found that ro¬ 
tational flows formed in each of the polarity concentrations soon 
after the intersecting flux tubes became vertical. Furthermore, 
the rotational movement of sunspots in a 3D MHD simulation 
has been examined by Fan (2009) where significant vortical mo¬ 
tions developed as a torsional Alfven wave propagated along 
the tube. [Fan[ ( [2009[ ) noted that the rotation of the two polari¬ 
ties twisted up the inner fieldlines of the emerged field, thereby 
transporting twist from the interior portion of the flux tube to 
the stretched coronal loop. In this simulation, a cylindrical flux 
tube is insert ed into the s olar interior using the cylindrical model 
developed by [Fan[ ( [2()01[ ). 

[Hood et aL] ([2009[ ) perform a complementary simulation to 
that of Fan ( [2001[ ), replacing the initial cylindrical flux tube with 
a toroidal tube. A common shortcoming of the cylindrical model 
in simulations without convective flows is that the axis of the 
tube never fully emerges. Altering the geometry of the flux tube 
to a curved shape allows for the axis of the tube to rise into the 
corona. Current theories suggest that the Sun’s magnetic field 
is created in the solar interior by dynamo action. Hence, a half¬ 
torus shaped flux tube imitates a field anchored deeper within the 
solar interior. The rotation of sunspots has not yet been investi¬ 
gated using this toroidal model; this is what we aim to study. 


In the present paper, we perform a resistive 3D MHD simu¬ 
lation of a toroidal flux tube placed in the solar interior and track 
its emergence through the photosphere and lower atmosphere. 
We study the role of rotational flows at the photosphere while 
also investigating the transport of magnetic helicity and energy. 
In addition, we explicitly calculate the angles of rotation in our 
experiment to directly compare with observations. Our main aim 
is to demonstrate that the interior magnetic field is untwisting as 
the tube emerges causing an injection of twist into the atmo¬ 
sphere as well as a rotation of the sunspots. This will be accom¬ 
plished by performing a detailed study investigating quantities 
relating to both the magnetic field and plasma. 


The remainder of the paper is structured as follows. In Sec¬ 
tion we describe the MHD equations used in our model as 
well as outlining the numerical approach employed. We also de¬ 
scribe the initial setup of our experiment, detailing both the ini¬ 
tial atmosphere and the model of the sub-photospheric magnetic 
field inserted in the solar interior. In Section the simulation 
results are presented, with emphasis on the rotational motions 
that develop within the two polarities on the photospheric plane. 
Our analysis also includes an investigation of the sunspot ro¬ 
tation angle, the flow vorticity at the photosphere, the current, 
relative magnetic helicity, magnetic energy and twist. Finally, in 
Section]^ we summarise the conclusions and discuss the impli- 
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cations of this simulation and future projects that stem from this 
work. 

2. Model setup 

In this section, we outline the numerical setup of our experiment, 
detailing the initial atmospheric and magnetic field configura¬ 
tion. 


in the x and y directions and from -25 to 75 in the z direction 
in dimensionless units. This corresponds to a physical size of 17 
Mm^. The boundary conditions are periodic on the side walls 
of the computational domain and the top and bottom boundaries 
are closed with v = 0. For all other variables we set the normal 
derivatives to zero. As a consequence, the magnetic field is fixed 
on the bottom boundary. 


2.1. Model equations and numerical approach 


For the experiment performed in this paper, we numerically 
solve the 3D time-dependent resistive MHD equations, as de¬ 
scribed below in non-dimensionalised form: 


— = -pV • v. 


Dp 

Dt 

Dv 

DB 


-(V X B) X B - -Vp ■ 
P P 


gz-\- 


V-T 


(B • V)v - B(V • v) - V X (pV X B), 

_^V-v + -f+ 

P P P ’ 


with specific energy density, 
P 


6 = 


(r- i)p’ 

and the electric current density, 

j = VxB. 


( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 


The basic dimensionless quantities used in these equations are 
the density p, time t, velocity v, magnetic field B, pressure p, 
gravity p, electric field E, specific energy density 6, resistivity rj 
and temperature T. The ratio of specific heats, y, is taken to be 
5/3. The viscosity tensor is denoted by T and the contribution 
of viscous heating to the energy equation is represented by Qvisc- 
The difference between simulations with and without a small vis¬ 
cosity are found to be negligible. The variables are made dimen¬ 
sionless against photospheric values. These values are: pressure, 
Pph = 1.4 X 10"^ Pa; density, Pph = 3 x 10“"^ kg m“^; temperature 
Tph = 5.6x 10^ K; pressure scale height, = 170 km; velocity, 
i;ph = 6.8 km s"^; time, ^ph = 25 s and magnetic field, 5ph = 1300 
Gauss. The dimensionless resistivity p has been taken as uniform 
set at a value of 0.005 in our experiment. All quantities from now 
on will be dimensionless, unless units are provided. In order to 
reach the true quantities with physical units, the dimensionless 
quantities should be multiplied by their photospheric values as 
given above. 

The code used to simulate the emergence process is a 3D La- 
grangian remap code, LareSd ( [Arber et al.|2001| ). The code uses 
a staggered grid and is second order accurate in both space and 
time. The FARE code can be divided into two main steps; the 
Lagrangian step and the remap step. In the Lagrangian step the 
equations are solved in a frame that moves with the fluid. This 
causes the grid to be distorted, so in order to put the variables 
back on to the original (Eulerian) grid, a remap step is used. 
The code accurately resolves shocks by using a combination of 
artificial viscosity and Van Leer flux limiters. This code also in¬ 
cludes a small shock viscosity to resolve shocks and the associ¬ 
ated shock heating term in the energy equation. 

The equations are solved on a uniform Cartesian grid (x, y, 
z) comprised of 512^ grid points. The box spans from -50 to 50 


2.2. Initial atmosphere 


The initial stratification of the atmosphere is the same as many 
previous flux emergence studies, for example |Hood et al!] ( |2009 ). 
The computational domain is split into four regions: the solar in¬ 
terior (z < 0); the photosphere/chromosphere 0 < z < 10; the 
transition region 10 < z < 20 and the lower corona z > 20. The 
solar surface is taken to be the plane at z = 0. The stratification is 
uniform across the horizontal plane and varies only with height z. 
The solar interior is taken to be marginally stable to convection 
by assuming constant entropy in this region. This assumption 
seems appropriate as the focus of this experiment is the evolution 
of the magnetic field. The photosphere/chromosphere is taken as 
an isothermal region with a dimensionless temperature of unity 
by design. The temperature distribution in the transition region is 
a power-law profile to model the steep temperature gradient ob¬ 
served here. Lastly, the isothermal corona is set at a temperature 
150 times that of the photosphere (Tcor = 150). To summarise, 
the non-dimensionalised temperature profile is specified as 



z < 0, 

1 

0 < z < 10, 

150(^) 

10 < z < 20, 

150 

z>20. 


The pressure and density are then calculated by numerically 
solving the dimensionless hydrostatic balance equation ^ = 
-pg and the dimensionless ideal gas law. The resulting loga¬ 
rithms of the initial temperature, density and pressure of the 
stratified background are shown in Figure 



Fig. 1: Initial stratification of model atmosphere. The initial pro¬ 
files are plotted on a log scale against height where red denotes 
the temperature distribution, black denotes the density and blue 
represents the magnetic pressure in the solar interior and the gas 
pressure throughout the domain. 
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2.3. Initial magnetic field 

We choose to leave the atmosphere unmagnetised by neglecting 
an ambient field and concentrating on the sub-photospheric field. 
We insert a magnetic fiux tube into the hydrostatic solar interior. 
As we want the setup to remain in equilibrium, and the tube to 
be in force balance, we require 

-Vp+jxB+pg = 0. 

Given that the background environment is in hydrostatic pressure 
balance, this reduces to 


^ Pqxc — j ^ C 7 ) 

where pexc is the pressure excess such that the gas pressure in the 
tube, pt, is defined to be + pexc where ph is the background gas 
pressure. Explicitly, we have split the gas pressure into a back¬ 
ground pressure component that balances gravity and a second 
component that balances the Lorentz force. The next obstacle is 
to choose the form of magnetic field to prescribe in the solar inte¬ 
rior. The interior magnetic field cannot be observed therefore we 
choose simple models for the sub-photospheric field to initiate 
emergence. We choose to place a toroidal fiux tube with twisted 
fieldlines in the solar interior. The toroidal model we utilise was 
derived fully in [Hood et al.| ( |2009] ). Compared with the standard 
choice of a cylindrical tube, this models the emergence of the top 
part of an Q-shaped loop that is rooted much deeper in the solar 
interior. Here, we simply outline the main steps of the derivation. 

First, we express our original cartesian coordinates (x, y, z) 
in terms of cylindrical coordinates (R, 0, x), given that y is the 
coordinate describing the direction of the axis of the tube, and z 
denotes the height from the solar surface as previously. We note 
that 

=y^ -^(z- Zbase)^, with y = Rcoscp and z - Zbase = R sin0, 


where Zbase is the value of z at the base of the computational 
domain. The magnetic field is then expressed, in terms of the 
fiux function A = A(R, x), as 


B = VA X V0 -h ^000 

IdA IdA 

where A is constant along magnetic fieldlines. In this derivation, 
we assume that the magnetic field is rotationally invariant, i.e. 
independent of 0 . We note that this form of magnetic field auto¬ 
matically satisfies the solenoidal constraint (V • B = 0). Inserting 
the magnetic field, B, in terms of the fiux function. A, into equa¬ 
tion 0, noting pexc = E(A(7?, x)) and RB^p = G(A(R, x)) along 
with some algebraic manipulation yields. 


1 

R 


Id^A 

r'^ 


1 db(h 


dp^ 

dA 


VA, 


where We note that all the terms are in the direction 

of VA, therefore the Grad-Shafranov equation is of the form 


lldA\ d^A ^db^ _ 

^dR (/? 5/?) 5x2 ^^<1’ dA 


2 dpQ^Q 

dA 


= 0 . 


( 8 ) 


Following the strategy used by |Hood et al. ( 2009| ), we now 
convert our system to a local toroidal coordinate system (r, 6, 0 ), 
such that 


(R- rI), with R- Rq = r cos 6 and x = -r sin 6, 
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where Rq is the major axis of the toroidal loop and a is the mi¬ 
nor radius of the toroidal loop. Equation 0 can be re-written 
in these local coordinates by assuming a « Rq and in turn 
r « Rq, i.e. assuming that the minor radius of the torus is much 
smaller than the major radius. We can expand the solution in 
powers of g/Rq such that to leading order, equation becomes 


Bed , 1 d{R^BllRl) dp,,, 

Jr 


= 0 . 


This has exactly the same form as the standard cylindrical equa¬ 
tion found in |Archontis et al.] ( |2004| ). We can therefore choose 
the solutions to be the same as the cylindrical flux tube. Specifi¬ 
cally, 


= and Bg = arB^ = a^Bore-d'^", 

K K 


where is the axial field strength and a is the twist. The local 
toroidal field can also be approximated to O(a/Ro). Again using 
R = Rq rcosO ~ 7?o + OialRo) this approximates the local 
toroidal field to 


^0 ~ Boe and Bq ~ aBore . 


With these approximations to the local toroidal field, the pres¬ 
sure can be calculated by exact comparison with the cylindrical 
model. 


Pexcir) = _ 2jr^ - 2). 

This balances equation 0 and forces the fiux tube to sit in equi¬ 
librium. However, to initiate the emergence we introduce a den¬ 
sity deficit as given by 


Pdef(0 ~ 


4T(z) 


_ 2 q ,2^ _ 2)_ 


using the temperature profile specified in Section |2.2| This 
makes the fiux tube lighter than its surroundings and allows it 
to begin rising. 

Eastly, we need to re-express the magnetic field in terms of 
cartesian coordinates as this is the form that we require for the 
input of our MHD simulations. Before we can do this directly, 
we must first express the field in cylindrical coordinates (/?, 0 , x) 
as 


Br = -Be{r) sin 6 = Be{r) -, 

r 

R — R 

Bx = -Bg(r) CO?, 0 =-Bg{r) - 

r 

Similarly, this cylindrical field can then be converted into carte¬ 
sian coordinates to set the magnetic field as 


Bx = -Bg{r )-— 

r 

By = 

B, = 
where 

^0 = BQe~^^^^^ and Be = arBcj, = aBQre~^^^^^. 

We note that this corrects an error in [Hood et ^ ( |2009| ) where 
the sign of Br and B^ were interchanged. Fortunately, this error 
was not significant as it is equivalent to using an initial twist, a, 
of the opposite sign. 
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2.4. Parameter choice 



Pseudocolor 
Var: Density 
-34.28 


I 


Pseudocolor 
Van Temperature 
-r 150.0 


28.23 

-5.313 

1.000 


Fig. 2: Summary of initial setup. The background density distri¬ 
bution is shown on the right wall, the temperature distribution on 
the back wall, a selection of fieldlines are shown in purple and 
an isosurface of the magnetic field (|B| = 1) is overplotted. 


In this experiment, we set the magnetic field strength at the 
apex of the tube as .So = 9 (11700 G) and the twist SiS a = 0.4 
(right-hand twisted). This means that each fieldline undergoes 
an angle of 0.4 rad per unit distance, and is equivalent to approx¬ 
imately three full turns of twist in the initial sub-photo spheric 
field. The base of the computational domain is set at z = -25. 
The major radius of the torus is So = 15 and the minor radius 
is a = 2.5. The total fiux threading a cross section of the tube 
is 6.6 X 10^^ Mx, typical of a small active region or a large 
ephemeral region. The initial set-up of the experiment is sum¬ 
marised in Figure 


In order to i nitiate this instab ility, a criterion must be satisfied 
as derived by | Acheson] ( [ 1 979| ) . Typically, the onset of this insta¬ 
bility occurs when the plasma p drops to one. At this stage, the 
magnetic pressure excee ds the gas pressure and t he fi eld expands 
into the atmosphere. See |Archontis et al.] ( |2004| ) and |Hood et al.| 
( 2012| ) for a full description of this instability. 



(b) t = 100 


Fig. 3: Visualisation of the field in the interior at times r = 40 and 
r = 100 respectively as traced from the lower negative footpoint 
(left). A movie of this figure is included in the electronic version. 


3. Analysis 

In order to examine and quantify the rotation of the sunspots in 
our simulation, we undertake a detailed investigation of a num¬ 
ber of quantities relating to both the plasma and magnetic field in 
an attempt to demonstrate the untwisting of the interior field and 
rotational flows at the photosphere. First of all, we briefiy de¬ 
scribe the overall evolution of the field detailing the main phases 
of emergence. 

3.1. Evolution of magnetic field 

The fiux tube begins to rise buoyantly to the photosphere due 
to the initial density deficit introduced and continues to rise be¬ 
cause of the decreasing temperature stratification in the interior. 
However, at the photosphere, the plasma is stably stratified with 
a constant temperature and the fiux tube is no longer buoyant. 
Therefore, the magnetic field finds another way to rise and ex¬ 
pand into the corona, namely the magnetic buoyancy instability. 


In Figurej^ we have included two figures illustrating the evo¬ 
lution of the magnetic fieldlines as the fiux tube emerges con¬ 
centrating on the interior portion of the field. We have traced 
three fieldlines from (0, -14, -25) (blue), (0, -15, -25) (black) 
and (0, -16, -25) (red) respectively. Initially the fiux tube has 3 
full turns of twist in total. At t = 40, the fiux tube reaches the 
photosphere and the legs start to straighten. At this time, there is 
an emerged section of fiux tube with about half a turn of twist 
that will subsequently expand into the corona. However, there is 
still a considerable amount of twist submerged, approximately 
a full twist in each leg. In the subsequent evolution there is a 
definitive unwinding of the submerged twist leading to a final 
state in which there is virtually no twist in the interior portion of 
the tube as evidenced by the visualisation of the field ait = 100.. 
We suggest that this may be governed by torsional Alfven waves. 
The release of sub-photospheric twist as a mechanism for the 
propagation of a torsional Alfven wave is investigated later. 
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3.2. Driver of rotational motion 


Before we examine various quantities at the photosphere in an 
effort to demonstrate the untwisting of the interior field, we 
discuss the underlying reason for a rotational movement at the 
photosphere when a twisted magnetic structure emerges. The 
basic cause for this rotational motion is the behav iour of the 
Lorentz force. Following the explanation given by [Cheung & 
jlsobej ( |2014| ), we consider a circular closed curve lying on the 
photospheric plane enclosing some point P denoting the loca¬ 
tion of the maximum of B^. We note that the torque is the ten¬ 
dency of a force to rotate an object about an axis, and is given by 
T = r X F where r is the displacement vector and F is the given 
force. If we consider the torque due to the various forces acting 
on the plasma and magnetic field through the surface confined by 
this closed contour, we find that the magnetic tension is the only 
force that provides a net torque. That is, the torque contributions 
from the magnetic pressure and gas pressure forces through this 
surface vanish and any non-zero torque is a result of the mag¬ 
netic tension. Explicitly, 


'^mp 


'^mt 


JJ^ r X V(-Pgas) • dS = 0, 




(B V)B dS, 


where r is the displacement vector of a point on the curve from 
P. This is in fact a general result that can be proved for any force 
of the form F = V/. To demonstrate this we have calculated the 
torque due to the magnetic tension and magnetic pressure within 
a circular contour of radius a (2.5) surrounding the location of 
the maximum of B^, as displayed in Figure In this case, it 
is clear that there is no contribution from the magnetic pressure 
force. Hence, we speculate that the driving motion of the rotation 
at the photosphere may be governed by the unbalanced torque 
produced by the magnetic tension force. This is characteristic of 
an Alfven wave. 



Fig. 4: Surface integral of torque due to the magnetic tension 
force (red) and the magnetic pressure force (blue) within a cir¬ 
cular contour of radius 2.5 around a point P corresponding to the 
maximum of the vertical magnetic field. 


3.3. Rotation angle 

As our previous analyses suggest that the interior portion of 
the flux tube is untwisting and significant rotational motions 


develop within the sunspots, we have again traced three field¬ 
lines from the base to the photosphere in an attempt to visu¬ 
alise this motion and quantify the amount by which the fieldlines 
have rotated. The axis of the fiux tube has been traced from the 
lower negative footpoint as well as two fieldlines either side of 
the axis in the ^-direction, i.e. we have traced fieldlines from 
(0, -14.5, -25) (blue), (0, -15, -25) (black) and (0, -15.5, -25) 
(red). A schematic of the traced fieldlines is shown in Figure [5a| 
and Figure for times ^ = 40 and times ^ = 80 respectively. 
These figures clearly demonstrate the rotation of the fieldlines 



(a) t = 40 (b) t = 80 


Fig. 5: Visualisation of the axis of the fiux tube (black fieldline) 
as well as two fieldlines (red and blue) spaced either side of the 
axis for comparison at selected times. A movie of this figure is 
included in the electronic version. 

threading the sunspots in a visual manner. Examining Figurej^ it 
appears that both the red and blue fieldlines have rotated through 
an angle of at least tt radians over 40 time steps. 

In order to follow the fieldlines undergoing this rotation, 
we have traced the v and y coordinates of the locations of the 
red, black and blue fieldlines as they pass through the photo- 
spheric plane. The trajectories of these fieldlines are shown in 
Figure Initially, the three fieldlines drift outwards in a line 
as the sunspots separate. Later, the motions slow down and ro¬ 
tation develops. This is not particularly clear given the transla¬ 
tional aspect of the motion as the sunspots separate. To remove 
this feature of the motion, we have subtracted off the location 
of the black axis from the position of the blue and red fieldlines 
in Figure This gives us the relative position of the blue and 
red fieldlines and indicates that the blue and red fieldlines rotate 
clockwise around the black axis by an angle of almost 360°. 


120 



with black axis location sub¬ 
tracted. 

(a) Fieldline trajectories. 


Fig. 6: The trajectories of the fieldlines as they pass through the 
photospheric plane coloured with increasing intensity as time 
progresses. The colour scale on the right shows the times dur¬ 
ing the evolution. 
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With the X and y coordinates of the intersections of selected 
fieldlines through the photosphere, we can calculate the angle of 
rotation using 

. , yo ~ ysLxis 

tan (j) = - 

Xq “ -^axis 

where Xaxis and ^axis are the x and y coordinates of the axis of 
the tube (black fieldline) and xq and yo are the coordinates of 
the fieldline we are investigating, i.e. the red or blue fieldline. 
In Figure [7] a schematic has been included to help us visualise 
the meaning of the angle 0. Through the use of equation the 
angle 0 has been calculated for both the red and blue fieldline as 
displayed in Figure As the two chosen fieldlines were initially 



(9) 


(xo, yo) 



Fig. 7: Representation of angle 0. 

equally spaced on either side of the axis, the rotation angles are tt 
out of phase at the beginning. Both fieldlines undergo a rotation 
of between InjA and In over 90 time steps. More precisely, the 
red fieldline undergoes a rotation of 340° and the blue fieldline 
undergoes a rotation of 353°. 



Time 


Fig. 8: Evolution of angle of rotation 0 for both the red and blue 
fieldlines as depicted in Figures [5a| and [5b| 

Next, we consider the temporal rate of change of the angle 
of rotation, i.e. d0/dr, as shown in Figure This illustrates that 
different regions of the sunspot are rotating at slightly different 
rates. There is an initial peak in the size of the rotation rate of 
both fieldlines of between -0.15 and -0.2. This peak in rotation 
rate occurs at about t = AA for the red fieldline and at about 
t = 62 for the blue fieldline. The rate of rotation diminishes as 
the experiment proceeds until it reaches zero indicating that the 
fieldlines have essentially stopped rotating. 

This rotation is investigated further by checking if the as¬ 
sumption of solid body rotation is reasonable, i.e. that the ro¬ 
tation angle does not depend on the radius of the fieldline. We 
have concluded from above that different areas of the sunspot 


Fig. 9: Evolution of rate of change of the angle of rotation 0 for 
both the red and blue fieldlines as shown in Figures [5a| and [5b| 



(a) Blue fieldline traced from (b) Red fieldline traced from 
(0, -14.5, -25). (0, -15.5, -25). 

d0 ojz 

Fig. 10: Comparison of terms — (solid line) and — (dashed 
line) for both the blue and red fieldlines respectively. 


are rotating at slightly different rates. If we assume the rotation 
is a solid body rotation then it follows that the velocity in the 0 
direction, at radius R, is given by 


„d<P 

="if- 


and the z-component of the vorticity is given by 




iA 

RM 


{Rv^) 



where we have assumed that 0 does not depend on R. We can 
therefore relate the vertical vorticity and the rate of change of 
the angle 0 by 


dt 2 ■ 


( 10 ) 


To check if our assumption is valid, we can investigate equa¬ 
tion ( p^ by plotting d0/dr and tu^/2 for both the red and blue 
fieldlines. In both panels of Figure [T^ the two terms balance 
each other well suggesting that equation ( p^ is valid and the ro¬ 
tation angle may not have a large dependence on the radius from 
the axis of the tube. 


3.4. Vorticity 

In another attempt to demonstrate the rotational movement at the 
photosphere, we analyse the plasma motion on the photospheric 
plane. Hence, the vorticity, calculated as the curl of the plasma 
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velocity, is examined as this quantifies the rotation of the plasma. 
As we are concerned with the rotation of sunspots, the vertical 
component of the vorticity is the quantity of interest, as it mea¬ 
sures the rotation in the x - ^ plane. This is expressed as. 


= (V X v)^ = 


dx 


dvx 
dy ■ 


A positive refers to a counter-clockwise motion while a nega¬ 
tive coz refers to a clockwise motion. In our experiment, the initial 
field is right-hand twisted; that is, if we were to create this field 
from a straight field, both footpoints would have been rotated 
in a clockwise motion. In other words, the fieldlines are wound 
counter-clockwise in the positive footpoint and clockwise in 
the negative footpoint. 

In order to visualise the field, let us consider the schematic of 
the fieldlines soon after the field has intersected the photosphere 
as shown in Figure At this point, the legs of the tube have 
started to straighten out and almost resemble a vertical cylindri¬ 
cal tube originating at z = -25 and intersecting the photosphere 
at z = 0. As the magnetic field is fixed at the base of the com¬ 
putational domain, a clockwise rotation at the photosphere is re¬ 
quired in order to unwind the interior portion of the field. The 
atmospheric field, on the other hand, is twisted up by a clock¬ 
wise rotational motion on the two footpoints at the photosphere. 



-20 -10 0 10 20 


(a) t = 40 


I 0.20 
0.13 
0 .. 
o.ca 

- 0.07 
- 0.13 
- 0.20 

-20 -10 0 10 20 
(b) t = 80 



Fig. 11: Coloured contours of accompanied by line contours 
of at z = 0, the base of the photosphere, for the specified 
times. A movie of this figure is included in the electronic version. 


the inner side of each of the sunspots again indicative of shear¬ 
ing motions. 

Now that we have established a clear rotational velocity at 
the photosphere, the evolution of the vertical vorticity at the pho¬ 
tosphere is expressed in a more quantifiable manner. Following 
the method used by |Fan|pQ09| ), we achieve this by plotting the 
time variation of the mean vertical vorticity, averaged over 
the area of the upper sunspot where B^ is greater than 3/4 of its 
peak value in Figure [T^ Explicitly, 




N \ 

yk, z = 0) 


( 11 ) 


where Xk and yk are the x and y coordinates of the region satis¬ 
fying Bz > ^msix(Bz) and N is the number of points that satisfy 
this criteria. This has been compared to the lower sunspot with 
no notable difference in result. 



Time 


Fig. 12: Evolution of the mean vertical vorticity {ojz) averaged 
over the area of the upper polarity concentration where Bz is 
above 75% of the max of on z = 0 as described in equa¬ 
tion ( pTj ). 


A series of coloured contour plots of the vertical vorticity are 
displayed in Figureat times ^ = 40 and t = 80. For visualisa¬ 
tion purposes, line contours of Bz have been over-plotted to show 
the location of the sunspots. At the centre of each of the polar¬ 
ities, there is a concentration of negative oJz corresponding to a 
clockwise rotation. This concentration of oJz appears when the 
legs of the tube straighten and builds with time until it peaks at 
about r = 50 before decaying as the simulation progresses. This 
hints that there is some bulk rotation of the sunspots, similar to 
the sunspot rotations observed by [Brown et al.| p003| ) and [Yanj 
et al.| ( |200^ . This result complements our previous investigation 
of the rotation angle where we found fieldlines undergoing sig¬ 
nificant rotations. Interestingly the same sign of vertical vorticity 
is present in both concentrations of opposite polarity. As dis¬ 
cussed earlier, this has important consequences for the evolution 
of the field in both the atmosphere and interior. We suggest that 
these sunspot rotations are due to the untwisting of the field in 
the interior injecting twist into the atmosphere. Another interest¬ 
ing feature of the contour plots is the red streak of negative vor¬ 
ticity between the sunspots. Elongated streaks of vorticity most 
likely correspond to shearing motions rather than rotation sug¬ 
gesting that these are due to shear fiows between the sunspots. 
In addition, there are blue tails of positive vorticity located on 


Considering the average vertical vorticity at the photosphere 
in Figure the vorticity is consistently negative suggesting 
that the dominant motion is a clockwise rotation. We find that a 
clockwise vortical motion appears in each polarity source when 
the field reaches the photosphere. The vortical motion quickly 
rises to a peak at roughly r = 50 soon after the emerged field be¬ 
comes vertical and the photospheric footprints have reached their 
maximum separation. At this time rapid rotation commences. 
The horizontal velocity at this time is approximately 0.5 in mag¬ 
nitude which corresponds to a physical velocity of 3.4 km s“^ 
Soon after, the vorticity steadily begins to decline. This signif¬ 
icant clockwise rotation twists up the emerged fieldlines in the 
atmosphere transporting twist from the tube’s interior portion to 
its stretched coronal portion. As noted earlier, we speculate that 
this transport of twist is due to some form of torsional Alfven 
wave. 

As discussed earlier, [Min & Chaej ( |2009| ) conjectured that 
the observed rotation of sunspots due to flux emergence may be 
an apparent effect due to a twisted field rising and each fieldline 
appearing in a different position at the photosphere. However, in 
our experiment we have established a clear rotation in the plasma 
velocity suggesting that this may not be the case. To estimate 
the contribution to the rotation by apparent effects, we quantify 
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the vertical advection of the flux tube by averaging the vertical 
velocity in a similar fashion to the way in which we averaged the 
vorticity. To achieve an upper bound for our estimate, we take the 
vertical speed of the tube to be the maximum of and assume 
that the vertical leg has a full turn of twist at = 40 when the held 
intersects. This is equivalent to the held being advected vertically 
by 2.4 units by the end of the experiment, corresponding to a 
34.6° apparent rotation. This is an over-estimate for the apparent 
rotation angle as we have assumed the maximum velocity for all 
time, and yet, this is still signiflcantly smaller than the calculated 
rotation angle. To conclude, this helps us disregard this theory 
and explain the rotation solely by a dynamical consequence of 
the emergence of flux. 

Though the vorticity plots do give us a useful insight into the 
rotational properties of the plasma at the photosphere, further 
examination is necessary in order to quantify the untwisting of 
the interior held, and the transport of twist into the atmosphere. 

3.5. Current density 

A useful quantity when estimating the twist of the magnetic held 
is the current density, speciflcally the z-component, as denoted 
by 


H\dx dy I' 

as this is related to how twisted the held is in the x - y plane. 
We consider coloured contours of at a height half way down 
the solar interior and at the base of the photosphere, as shown 
in Figure From the coloured contours of in Figure it 



(c) z = 0 and t = 40 


(d) z = 0 and t = 80 


Fig. 13: Coloured contours of at the plane in the middle of 
the interior (z = -12.5) in the top panel and at the solar surface 
(z = 0) in the bottom panel for the specifled times, as well as line 
contours of for comparison of the size of sunspots. 


is clear that although the majority of each sunspot is positive or 
negative, the periphery of each sunspot is dominated by the op¬ 
posite sign of jz. As we insist that our initial sub-photo spheric 
flux tube is isolated and therefore surrounded by unmagnetised 
plasma, Faraday’s law demands that the flux tube must carry no 


net current. As the flux tube carries current inside due to its twist, 
a reverse current surrounds the sunspot to ensure a zero net cur¬ 
rent in this region. Focussing on the top panel of Figure[^ there 
is some evidence that the two concentrations of strong centred 
on the sunspots are depleting with time. Similarly, considering 
the bottom panel of Figure the concentrations of at the 
photosphere intensify when the held first emerges then diminish 
as the experiment proceeds. As signifies the twist of the held 
in the V - ^ direction and is related to the azimuthal magnetic 
held, a decrease in could indicate a decline in the amount of 
twist. The mechanism responsible for this decrease in requires 
further investigation. 



Fig. 14: Plot of the maximum value of against time for varying 
heights depicted in the key. 

The temporal variation of the maximum value of for spe¬ 
cific heights below the photosphere is displayed in Figure 
There is an initial increase in the maximum of for all heights 
due to the emergence of the held before a steady decline as the 
experiment proceeds. There are two possible explanations for 
this steady decrease in the vertical current. This could be caused 
by the expansion and stretching of the held as a result of emer¬ 
gence or by a decline in the amount of twist stored in the held. 
The stretching of the held results in a decline in the gradients of 
Bx and By and lowers j^. To evaluate the extent of the expansion 
of the held, we estimate the diameter of a contour of B^ as shown 
in Figure This plot shows the maximum y-separation of the 
contour of B^ = 1. Initially, the separation increases for heights 
near the photosphere as the flux tube buoyantly rises and the legs 
of the tube straighten. Later, there is very little change in the sep¬ 
aration of the legs for heights deep in the solar interior. There is 
however some expansion of the magnetic held for the photo- 
spheric height z = 0 and 5 units below the boundary as expected 
as the magnetic held expands into the low density atmosphere. 
This helps us disregard this cause and explain the decrease in 
in the solar interior solely by an untwisting of the held. Specif¬ 
ically, a decrease in implies a decline in the azimuthal held 
which corresponds to the interior portion of the held untwisting. 
However, a decrease in at the photosphere may be explained 
by the expansion of the held in this region. 

Finally, to analyse the current further, an examination of 
the total jz in each sunspot is necessary. As we noted from the 
coloured contours, although the centre of each sunspot is dom¬ 
inated by one sign of current, the outer boundary consists of 
reverse current. Hence, we cannot estimate the current in each 
sunspot by measuring the total positive or negative current. In- 
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Fig. 15: Line plot of the diameter of the = 1.0 contour as a 
function of time for varying heights depicted in the key. 


the mean vertical current generally declines in a linear man¬ 
ner. After an initial drop in ( 7 ^), there is a small increase as the 
field straightens out and the negative outer-boundary plays a less 
dominant role. Overall, there is a general decrease in (j^) sug¬ 
gesting a drop in the twist stored in the interior portion of the 
field as the field untwists. 

3.6. Magnetic helicity 

In order to quantify the transport of twist from the solar interior 
to the solar atmosphere due to the emergence and untwisting of 
the field in the interior, we have calculated the evolution of the 
relative magnetic helicity both above and below the photosphere. 
The magnetic helicity is a topological quantity describing how 
much a magnetic structure is twisted, sheared or braided. 

The relative magnetic helicity (to the refer ence field Bp) of 
the field B in a given volume V is given by ( [Berger & Field 
( |1984| ), [Finn & Antonsen| ( p^85] )) 

Hr = £(A + Ap)-(B-Bp)dV, (12) 



(a)z = -12.5 



Time 

(b)z = 0 

Fig. 16: Time evolution of the vertical current {j^) averaged over 
the area of the positive polarity fiux source where is greater 
than 75% of its maximum. 


stead, we estimate the current in the centre of the upper sunspot 
by averaging the vertical current over the area where the vertical 
magnetic field is greater than 3/4 of its maximum in a similar 
fashion to our calculation of the average vorticity. The temporal 
evolution of this quantity is depicted in Figure for the interior 
plane (z = -12.5) and the photospheric plane respectively. 

At the photosphere, the mean vertical current increases as 
the magnetic field emerges then steadily declines as the field ex¬ 
pands into the atmosphere. Lower down in the solar interior. 


where A is the vector potential of B (B = V x A), B^ is the 
reference potential field with the same normal fiux distribution 
as B on all bounding surfaces and A^ is the vector potential of 
Bp (Bp = V X Ap). The relative magnetic helicity is favoured 
over the magnetic helicity as it has been shown to be gauge- 
independent with respect to the gauges for A and Ap. This is 
a necessary condition for a physically meaningful definition of 
helicity. 

In order to calculate the relative magnetic helicity numeri¬ 
cally, we have tested and compared two app roaches; the ca lcu- 
lation of A and Ap using DeVore’s method (|DeVore||2000j) and 


the numerical procedure from Moraitis et al.|(|2014j). As the re 


suits were comparatively similar, we will discuss the calculation 
of the relative magnetic helicity from [Moraitis et al.] ( |2Q14| ). 

In calculating the potential field within the volume V = 
[x i,V 2 ] X [^i,j/ 2 ] X [zi,Z 2 ], the numerical procedure utilised 
in [Moraitis et al.[ ( [2014[ ) takes into account all boundaries within 
the finite volume. This is advantageous over DeVore’s method 
which is only valid for a semi-infinite space above a lower 
boundary. The potential field satisfies jp = V x Bp = 0 within 
y, thus implying Bp = -V 0 where 0 is a scalar potential. The 
solenoidal constraint V • Bp = 0 then implies that the scalar po¬ 
tential is a solution of Laplace’s equation = 0 in V. The con¬ 

dition that B and Bp have the same normal components along the 
boundaries of the volume translates to Neumann boundary con¬ 
ditions for 0 , i.e. d(pldh\Qy = - • B|^y. The Laplace equation is 

then solved numerically using a standard FORTRAN routine in¬ 
cluded in the FISHPACK library ( [Swarztrauber & Sweet|1979[ ). 

The original and potential fields are now stored for the given 
time step and desired volume. The next step is to calculate 
th e corresponding ve ctor potentials given the method proposed 
by [Valori et al.[ ( [2Ql^ . As the relative magnetic helicity is gauge- 
independent, we are free to choose the gauge A •£ = 0 throughout 
V so that the x and y components of B = V x A are integrated 
over the interval (zi , z) to 


f 


A = Ao-zx I B(x,y,z)dz, 


where Aq = A(x,y,z = Zi) = (Aox,Aoy,0) is a solution to the 
z-component ofB = V xA,i.e. 


SAqx 


dx 


dy 


= B^(x,y,z = Zi). 
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Following Valori’s method we choose the simplest solution to 
the above equation, given by 

1 n 

Aox(x,y,z = zi) = -- I B^(x,y',z = z\)dy', 

1 

Ao^(v,i/,z = Zi) = -J = z\)&x'. 

Similarly, the vector potential of the potential field is calculated 
using 

Ap=Ao-zx r 'Rp{x,y,z')dz!, 

where we have noted that A^o = Aq as B and share the same 
normal component on the boundary at z = zi. 

An alternative solution for the vector potentials can be ob¬ 
tained if we use the top boundary, i.e. integrating over the inter¬ 
val (z, zi) as 


rz2 

A = Ao+zx I B(x, ^,z') dz'. 

This has been checked for comparison and there is no notable 
difference between the two solutions. 



Time 


Fig. 17: Evolution of the relative helicity Hr when calculated 
above z = 0 in the solar atmosphere. 


With the corresponding potential and vector potentials for 
the magnetic field, we can calculate the relative magnetic helic¬ 
ity within different sub volumes of the total simulation domain. 
First, we calculate the atmospheric magnetic helicity above the 
photosphere, as shown in Figure [T^ As expected, the atmo¬ 
spheric helicity is zero until ^ = 20 when the magnetic field first 
emerges. Subsequently, there is a linear increase in helicity as 
twist is steadily injected into the atmosphere by the emergence 
of flux and rotational motions at the photosphere. By the end 
of the experiment, the magnetic helicity injected into the atmo¬ 
sphere has reached 3.6 x 10^^ Wb^ (3.6 x 10^^ Mx^), typica l of a 
small event. Th is can be compared with observations where [Mm| 


& Chae 


(2009) quote a helicity transport of 4 x 10^^ Mx^ when 


considering a much larger active region. To investigate the nor¬ 
malised magnetic helicity that is transported to the atmosphere, 
we divide by and find that the evolution follows the same 
linear monotonic shape and reaches a value of 0.83 hy t = 120 
corresponding to almost one full twist of the flux in the original 
tube. 

To further investigate the magnetic helicity in the atmosphere 
we have calculated the time derivative of the magnetic helic¬ 
ity above the photosphere, both numerically using finite differ¬ 
encing and analytically. This helps us to understand the main 



Fig. 18: Evolution of the time derivative of the relative helicity 
Hr when calculated above z = 0 in the solar atmosphere us¬ 
ing equation ( p^ . The total time derivative (black solid line) is 
split into the dissipation term (purple solid line), the surface cor¬ 
rection term (yellow solid line), the shear term (red solid line) 
and emergence term (blue solid line). The dashed black line is 
the derivative of the curve from Figure [^calculated using finite 
differencing. 


sources contributing to the production and depletion of helicity. 
Magnetic helicity is mainly contributed to by vertical flows that 
advect twisted magnetic flux into the corona and by surface flows 
that shear and twist magnetic fields (|Berger & Field|1984|). The 
rate of change of relative magnetic helicity. Hr, can be evaluated 
analytically by differentiating the expression in ^V2) (|Berger & 
|Field|1984) , 

dHr 
dt 


-ly J'j BdV + 2rjJ' ((Ap x j) • n d5 
2 J [(Ap • v)(B • n) - ((Ap • B)(v • n)] dS. 


(13) 


The first term on the right-hand side of equation ( [T^ relates 
to the depletion of helicity by internal dissipation (dissipation 
term), the second corresponds to a surface correction to the re¬ 
sistive dissipation (surface correction term) , the third relates to 
the generation of helicity by horizontal motions of the boundary 
(shear term) and the last corresponds to the injection of helicity 
by direct emergence (emergence term). Let us consider the rate 
of change of atmospheric helicity in Figure From the point 
the field emerges until r = 45, the helicity flux due to emergence 
(blue solid line) dominates the change in helicity. Later, the hor¬ 
izontal shearing and rotational motions at the photospheric foot- 
points (red solid line) are the primary sources of helicity change, 
in agreement with ( [Fan 2009| ). The change in helicity due to 
internal helicity dissipation (purple) and the surface correction 
term (yellow) are much less significant and do not contribute to 
the overall change in helicity in the atmosphere. The derivative 
of Hr from Figure [17] is over plotted for comparison. The two 
approaches agree very well indicating that numerical effects are 
kept to a minimum. Care must be taken when considering the 
rate of change of helicity given in equation ( [T3] ). See |Pariat et al.| 
( |2015| ) for the full derivative including additional terms. Clearly, 
from Figure[^ the additional terms are not important in this par¬ 
ticular case as the flux throug h the surfaces clos ely follows the 
time derivative of the helicity. |Pariat et ar] ( |2Q15| ) also notes that 
precaution must also be taken when dividing the helicity flux into 
individual terms as although their sum is gauge-independent the 
individual terms are not, hence, limiting their physical meaning. 

Now that we have established a clear increase in the helicity 
in the atmosphere, we analyse the helicity in the solar interior. 
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As the toroidal flux tube has shown clear signs of untwisting, 
we expect that the magnetic helicity in the interior will decrease 
due to both the emergence of magnetic flux and the rotational 
movements on the photospheric boundary. The calculation of 
the relative magnetic helicity and the change in helicity below 
the photosphere are shown in Figure [T^ As expected, the inte- 



(a)//. 




(a) Ffj-ee 



(b)Fp 

Fig. 20: Evolution of (a) the free energy when calculated above 
z = 0 and (b) the vertical Poynting flux through the surface z = 0. 
The total Poynting flux is split into the emergence term (blue), 
the shear term (red) and the resistive term (purple) as deflned in 
equation ( p^ . 


(b) dHJdt 

Fig. 19: Evolution of the relative helicity and rate of change of 
helicity when calculated below z = 0 in the solar interior. 


rior relative magnetic helicity decreases throughout the experi¬ 
ment. Initially, the decrease in magnetic helicity is dominated by 
internal helicity dissipation. Later, at ^ = 20, the held emerges 
from the interior through to the atmosphere resulting in a sharper 
decrease due to both the emergence of magnetic flux and the ro¬ 
tational movements on the photosphere extracting helicity from 
the interior as the flux tube untwists. 


3.7. Magnetic energy 


Now that we have demonstrated the transport of magnetic helic¬ 
ity from the interior portion of the domain to the coronal portion, 
we consider the magnetic energy and its distribution across the 
domain. In order to understand the amount of magnetic energy 
available for solar eruptive events, we consider the free magnetic 
energy relative to the potential held. That is, we calculate the 
excess magnetic energy after subtracting off the energy stored 


in the potential held, i.e. £’free 


J B^/2dV-J By 2 dV. The 


evoluti on o f the free magnetic energy above z = 0 is shown in 
Figure [20a| The excess energy above z = 0 builds from the time 
the held first emerges. To investigate the contributions to mag¬ 
netic energy by flux through the photospheric boundary we con¬ 


sider the Poynting flux through z = 0 as given by, 

Fp= I B^Vzdxdy- I (y'B)Bz dxdy-vrj I (jxB)kdxdy. 

Jz=0 *Jz=0 Jz=0 

(14) 

The first term contributing to the Poynting flux in equation ( p^ 
corresponds to the contribution by vertical flows owing to emer¬ 
gence, the second denotes the generation of magnetic energy by 
shearing/rotational flows and the third term is a result of resistive 
effects. The rate of increase of energy is largest during the initial 
stages of emergence due primarily to the emergence term. How¬ 
ever, later the shear term (attributed to by rotational motions) is 
the main contributor to magnetic energy increase in the atmo¬ 
sphere. This pattern corroborates the trend that appeared in the 
helicity flux whereby vertical flows dominate the flux initially 
and horizontal flows become important later. Keeping in mind 
the precaution in Section |3.6| about the individual terms of he¬ 
licity flux, the behaviour of the Poynting flux helps us to believe 
that the helicity flux trend may have physical meaning. At the 
end of the experiment, the free magnetic energy transported to 
the atmosphere has reached 8.2 x 10^^ J (8.2 x 10^^ ergs). 


3.8. Force-free parameter 

Now that we have considered the behaviour of the plasma flows, 
current and magnetic helicity and energy, a proxy for the local 
twist is presented. Consider the quantity, a, normally referred 
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Tlme=40.0018 


(a) t = 40. 



Time=100 


(b) t = 100. 



In order to investigate this expression, we trace a along field¬ 
lines as a tool to help us visualise the distribution of twist along 
the field as shown in Figure Initially the field is coloured red 
to signify the field is highly twisted with a value of a greater than 
0.2. Rapid emergence results in a coronal magnetic fiux that is 
initially quite untwisted ( [Longcope & Welsch|2000| ) as shown by 
the low coronal a. A clear gradient in a develops along the field¬ 
lines from a high magnitude a in the interior to a lower magni¬ 
tude of a in the atmosphere as displayed in Figure [^|Fan|p00^ 
suggests that this gradient produces a torque that drives the ro- 
tational motions observed. This result has been proven by [Long- 1 
cope & Klapper| \\991) and it is suggested that this motion will 


continue until the gradient in a is removed ( [Longcope & Welsch| 
|2000| ). This agrees with our results. 

Later, after the rotational motions have set in at the photo¬ 
sphere, the interior field is left coloured blue with weak twist 
{a < 0.2) and the atmospheric field is threaded with a mixture 
of blue, red and white fieldlines. Thus, we can conclude that the 
field in the legs of the tube undergo a definitive untwisting mo¬ 
tion as the field is transformed from highly twisted to weakly 
twisted by the end of the simulation. Although a large portion of 
the atmospheric field is coloured blue due to the expansion and 
stretching of the field, some fieldlines are coloured red indicating 
that some highly twisted field has been t ransp orted into the at¬ 
mosphere. This is demonstrated in Figure [2^ where the twisted 
structure of the atmospheric field is shown. As a ~ 1 /L, an in¬ 
crease in the length scale results in a decrease in a. The length 
scales in the corona are much larger so a decrease in a in the 
atmosphere can be explained by the expansion and stretching of 
magnetic fieldlines. The length scales in the interior, on the other 
hand, are much smaller and so we explain a decrease in a by an 
untwisting of the field. 


3.9. Propagation of torsional Alfven wave 

As touched on earlier, our results indicate that the rotational mo¬ 
tions we observe may be governed by some form of torsional 
Alfven wave propagating downwards as a consequence of the 
transport of twist from the interior to the corona. The travel time 
for an Alfven wave to propagate from the photosphere to the base 
of the domain is approximately 20 time steps. This suggests that 
an Alfven wave would take approximately 40 time steps to travel 
down to the base, refiect and return to the photospheric plane. 
We propose that this downward propagating Alfven wave, first 
proposed by [Longcope & Welsch| ( |2000| ), untwists the magnetic 
field in the interior. The rotation will only slow down once the re- 
fiected wave has returned to the photosphere. This appears to be 
in fairly good agreement with Figure [T^ where the rapid rotation 
and large occurs from about r = 50 to r = 90. 

4. Conclusions/Discussions 


(c) t = 100. 

Fig. 21: Visualisation of magnetic fieldlines traced from both 
footpoints coloured by the parameter a such that red represents 
a strong twist (0.2 < a < 0.4) and blue denotes a weaker twist 
(0 < a < 0.2). 


to as the force-free parameter or sometimes the fieldline torsion 
parameter, 

O' = (V X B) • B/bI 


We have presented a 3D MHD numerical experiment of the 
emergence of a toroidal fiux tube from the solar interior through 
the photosphere and into the solar corona. Based on our detailed 
analysis, there is evidence that the interior magnetic field un¬ 
twists and the photospheric footprints undergo a rotation. This 
rotational motion acts to untwist the interior field fixed at the 
base and injects twist into the emerged atmospheric portion. Our 
analysis of the plasma vorticity at the photospheric plane re¬ 
veals that significant vortical motions develop in the centre of 
the sunspots. A definitive rotation of the sunspots is also demon¬ 
strated by tracking the fieldlines and calculating the rotation rate 
of the fieldlines threading the sunspot. Rotations of the order of 
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one full rotation (360°) are observed in our experiment. This is 
similar in magnitude to the angles of rotation reported in studies 
that concluded a direct relationship between swift sunspot ro¬ 
tation and en han ced eruptive activ ity ( [Brown et al.| ( |2003| ), [Yan 


& Qu (2007) and Yan et al.] ( 2009| ) etc.). However, the sunspots 


in our experiment rotate by one full rotation over the course of 
about 40 minutes, whereas [Brown et ah] ( [2003 j ) found sunspots to 
rotate through angles of the order of 200° over a period of days. 
The timescales in this experiment are clearly not in line with ob¬ 
servations. However, this is related to the size of our emerging 
active region (~ lOMm). If we scale up our experiment to a more 
typical active region, we expect the timescales to be much larger 
and in line with observations. This requires further investigation. 

The direct emergence of flux paired with the continual rota¬ 
tional motions at the photosphere transport magnetic energy and 
helicity from the solar interior to the corona. The magnetic en¬ 
ergy in the atmosphere reaches a value of 8.2x 10^^ J by the latter 
stage of the experiment. The Poynting flux of energy is split into 
contributions from horizontal shearing and vertical emergence 
terms. The initial flux of energy across the photospheric bound¬ 
ary is dominated by emergence but latterly the primary source is 
the horizontal shear. The rate of change of relative magnetic he¬ 
licity in the solar atmosphere also has two predominant sources; 
namely helicity flux due to emergence and helicity flux by ro¬ 
tational motions. In a similar trend to the energy, initially the 
predominant source of helicity is the emergence of the magnetic 
flux tube but later this is replaced by the flux of helicity due to ro¬ 
tational motions at the photospheric level. The magnetic helicity 
transported to the atmosphere reaches a value of 3.6 x lO^^Mx^. 
As well as the production of helicity in the atmosphere, we And 
a clear decrease in the magnetic helicity in the interior, support¬ 
ing our understanding that this portion of the held undergoes an 
untwisting motion as also evidenced by a clear decrease in in 
the interior. 

The appearance of rotational motions centred on both 
sunspots has been found before in other MHD simulations in¬ 
cluding [Magara[ ( [2006[ ) and [Fan[ p009[ ). [Fan[ ( [2009[ ) also ex¬ 
plained these rotations as a consequence of torsional Alfven 
wave propagation and established an increase in helicity in the 
atmosphere. Our work, however, explicitly discusses the effect 
that this rotational motion has on the interior portion of the held 
by establishing a depletion in the magnetic helicity stored in the 
interior segment of the domain and a drop in the vertical cur¬ 
rent in this region. We also show that the magnetic tension force 
may govern this rotational motion as it appears to produce an 
unbalanced torque that drives the rotation. Furthermore, we ex¬ 
plicitly rule out that the rotation observed may be an apparent 
effect, helping us to explain the rotation primarily as a dynam¬ 
ical result of the emergence of magnetic flux. In addition, we 
trace fleldlines from the base of the domain as they pass through 
the photosphere and explicitly calculate their angles of rotation 
which appear to be approximately in line with the angles cal¬ 
culated in observations. By considering the trajectories of these 
selected fleldlines, we And a remarkable visual representation of 
two fleldlines rotating in a circular motion around the axis (the 


centre of the sunspot), as shown in Figure 6b 


The simple hypothesis from [Longcope & Welsch[ ( [2000[ ) has 
been confirmed in this 3D resistive MHD model, namely that 
only a fraction of the current carried by a twisted flux tube 
will pass into the corona and that the propagation of a torsional 
Alfven wave at the time of emergence will transport twist from 
the highly twisted interior to the stretched coronal loops. The 
rotational motions we observe are a manifestation of the propa¬ 
gation of these waves. 


In a future paper we plan to perform a parameter study where 
we investigate the effects of the held strength and twist a on 
our analysis of the rotation of sunspots. This will allow us to de¬ 
termine if there is a relationship between the strength or twist of 
the sub-photospheric flux tube and the amount of vorticity in the 
sunspots or the rate at which the sunspots rotate. Another possi¬ 
ble avenue for future research is to add in an ambient magnetic 
held to the coronal portion of the simulation domain. This would 
add in a further realism and allow us to understand the effect of 
an ambient held to the rate of rotation, as well as understanding 
whether rotation can lead to an eruption in our experiments. 

Furthermore, we would like to try and understand why the 
rotation rate we calculate in our simulation is much larger than 
those calculated in observations. There are many possibilities for 
this discrepancy, including the size of active region which was 
discussed earlier. In addition, varying the strength or twist of the 
tube may change the time it takes for the flux tube to rise to the 
photosphere and he nce govern the rotation rate of the tube. The 
model presented by [Longcope & Welsch[ ( [2000| predicts that the 
level of rotation will depend on the rapidity of flux emergence so 
we plan to investigate how this affects the rotation. The length of 
time for the rotation may also be related to the depth at which 
the flux tube is anchored; this is another approach that requires 
investigation. 
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